home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QREMAIN.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  2KB  |  112 lines
  1. /*
  2. This is borrowed from ieee.c, where it is used for radix conversion.
  3. Hope the program didn't break in translation.
  4. It would execute very slowly if numerator >> denominator.
  5. Detecting whether n * a is exact could speed it up.
  6. */
  7.  
  8. /* Floating point remainder.
  9.    c = remainder after dividing b by a.
  10.    If n = integer part of b/a, rounded toward zero,
  11.    then qremain(a,b,c) gives c = b - n * a.  */
  12.  
  13. #include "qhead.h"
  14. #include "mconf.h"
  15. extern QELT qzero[];
  16. int normlz(), cmpm(), subm(), shdn1(), shup1(), shup16(), shup8(), mdnorm();
  17. int qcmp(), qclear(), qmov();
  18.  
  19. static QELT quot[NQ+1];
  20.  
  21. int qremain( a, b, c )
  22. QELT a[], b[], c[];
  23. {
  24. QELT den[NQ+1], num[NQ+1];
  25. QELT j;
  26. int i;
  27.  
  28. if( qcmp(a,qzero) == 0 )
  29.     {
  30.     mtherr( "qremain", SING );
  31.     qclear( c );
  32.     return 0;
  33.     }
  34. den[NQ] = 0;
  35. qmov( a, den );
  36. den[0] = 0;
  37. num[NQ] = 0;
  38. qmov( b, num );
  39. num[0] = 0;
  40.  
  41. /* Execute divide steps until num < den.
  42.    Least significant integer quotient bits left in quot[].  */
  43.  
  44. quot[NQ] = 0;
  45. qclear( quot );
  46. while( qcmp(num,den) >= 0 )
  47.     {
  48.     if( cmpm(den,num) <= 0 )
  49.         {
  50.         subm(den, num);
  51.         j = 1;
  52.         }
  53.     else
  54.         {
  55.         j = 0;
  56.         }
  57.     shup1(quot);
  58.     quot[NQ] |= j;
  59.     shup1(num);
  60.     num[1] -= 1;
  61.     }
  62. /* Normalize.  */
  63. while( num[2] != 0 )
  64.   {
  65.     shdn1( num );
  66.     num[1] += 1;
  67.   }
  68. i = 0;
  69. do
  70.   {
  71.     if( num[3] & SIGNBIT )
  72.       goto normok;
  73. #if WORDSIZE == 32
  74.     else if( (num[3] & 0xffff0000) == 0 )
  75. #else
  76.     else if( num[3] == 0 )
  77. #endif
  78.       {
  79.     shup16(num);
  80.     num[1] -= 16;
  81.     i += 16;
  82.       }
  83. #if WORDSIZE == 32
  84.     else if( (num[3] & 0xff000000) == 0 )
  85. #else
  86.     else if( num[3] & 0xff00 == 0 )
  87. #endif
  88.       {
  89.     shup8(num);
  90.     num[1] -= 8;
  91.     i += 8;
  92.       }
  93.     else
  94.       {
  95.     shup1(num);
  96.     num[1] -= 1;
  97.     i += 1;
  98.       }
  99.   }
  100. while( i<=NBITS );
  101. qclear(num);
  102.  
  103. normok:
  104.  
  105. /* Sign of remainder = sign of dividend. */
  106.  
  107. *num = *b;    /* Set sign of Remainder */
  108.  
  109. qmov( num, c );
  110. return 0;
  111. }
  112.